#ifndef _CGCM_PARAMS_SRC_
#define _CGCM_PARAMS_SRC_

#include "model.h"
#include "mexAux.h"
#include "netcdfaux.h"
#include <stdarg.h>
#include "readers.h"
#include <limits>
#ifndef SERIALQGCM
#include <omp.h>
#include "qgcmio.h"
#endif

class model_class::local { // all static functions
  public:
    static void profcheck_assert(bool check, const char * errormsg, ...)
    {
      if(!check){
        va_list(arg);
        va_start(arg, errormsg);
        vfprintf(stderr, errormsg, arg);
        va_end(arg);
        fprintf(stderr, "\n");
        throw cgcmexception(string("Profile/Parameter error."), -10);
      }
    };
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -//
    static void checkparams( model_class & mo ){
      const unsigned int & numfields = mo.get_global_param_by_name(gdb_str_NUMFIELDS).getiValue();
      profcheck_assert(
          numfields >= 0 && numfields <= 8, "Parameter '%s' must lie between %i and %i.", gdb_str_NUMFIELDS, 1, 8 );
      const REAL & threshold = mo.get_global_param_by_name(gdb_str_PQFINEDOMAINACCURACYTHRESHOLD).getfValue();
      profcheck_assert( threshold >= 0 && threshold <= 1.0, "Parameter '%s' must lie between 0 and 1.",
          gdb_str_PQFINEDOMAINACCURACYTHRESHOLD);

      const REAL & timestep = mo.get_global_param_by_name(gdb_str_TIMESTEP).getfValue();
      profcheck_assert(timestep > 0, "Parameter '%s' must be greater than 0.", gdb_str_TIMESTEP);
      const int & init_depth = mo.get_global_param_by_name(gdb_str_INITIAL_NESTING_DEPTH).getiValue();
      profcheck_assert(init_depth > 0, "Parameter '%s' must be greater than 0.", gdb_str_INITIAL_NESTING_DEPTH);

      const int & atmos_layer_1 = mo.get_global_param_by_name(gdb_str_ATMOSPHERE_LAYER_1_FIELD_NUMBER).getiValue();
      const int & atmos_layer_2 = mo.get_global_param_by_name(gdb_str_ATMOSPHERE_LAYER_2_FIELD_NUMBER).getiValue();
      const int & ocean_layer_1 = mo.get_global_param_by_name(gdb_str_OCEAN_LAYER_1_FIELD_NUMBER).getiValue();
      const int & ocean_layer_2 = mo.get_global_param_by_name(gdb_str_OCEAN_LAYER_2_FIELD_NUMBER).getiValue();
      const int & mixed_layer = mo.get_global_param_by_name(gdb_str_MIXED_LAYER_FIELD_NUMBER).getiValue();
      const int & topo_layer = mo.get_global_param_by_name(gdb_str_TOPO_LAYER_FIELD_NUMBER).getiValue();

      const int index_list [6] = { atmos_layer_1, atmos_layer_2, ocean_layer_1, ocean_layer_2, mixed_layer, topo_layer };
      const char * str_list[6] = { 
        gdb_str_ATMOSPHERE_LAYER_1_FIELD_NUMBER, gdb_str_ATMOSPHERE_LAYER_2_FIELD_NUMBER,
        gdb_str_OCEAN_LAYER_1_FIELD_NUMBER, gdb_str_OCEAN_LAYER_2_FIELD_NUMBER,
        gdb_str_MIXED_LAYER_FIELD_NUMBER, 
        gdb_str_TOPO_LAYER_FIELD_NUMBER 
      };
      for(size_t iter = 0; iter < 6; iter++){
        profcheck_assert( index_list[iter] > 0, "Parameter '%s' must be an integer larger than 0.", str_list[iter]);
        for(size_t u = 0; u < 6; u++){
          if(u != iter ){
            profcheck_assert( index_list[iter] != index_list[u], 
                "Parameter '%s' and parameter '%s' cannot point to the same field index.", str_list[iter], str_list[u]);
          }
        }
      }
    };
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -/
    static void mapdkprofile( model_class & mo){
      const vector<REAL> & kdomain = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN).getData();
      if( ! mo.does_prof_exist(gdb_prof_str_DK_DOMAIN) ){ // map it
        vector<REAL> dk;
        for(size_t i = 1; i < kdomain.size(); i++){
          dk.push_back( kdomain[i] - kdomain[i-1] );
        }
        mo.addprofile( new parameter(gdb_prof_str_DK_DOMAIN, dk) );
      }
    };
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -/
    static void mapkspacingprofile( model_class & mo ){
      const vector<REAL> & kc = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_CENTERED).getData();
      if( ! mo.does_prof_exist(gdb_prof_str_K_DOMAIN_SPACING) ){ // map it
        vector<REAL> kspacing;
        for(size_t i = 0; i < kc.size(); i++){
          kspacing.push_back( gdb_enum_UNIFORMSPACING ); // not using a string profile here because they are not supported
        }
        mo.addprofile( new parameter(gdb_prof_str_K_DOMAIN_SPACING, kspacing) );
      }
    };
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -/
    static void mapkcenteredprofile( model_class & mo){
      const vector<REAL> & kdomain = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN).getData();
      if( ! mo.does_prof_exist(gdb_prof_str_K_DOMAIN_CENTERED) ){ // map it
        vector<REAL> kc;
        for(unsigned int i = 1; i < kdomain.size(); i++){
          kc.push_back( 0.5*(kdomain[i] + kdomain[i-1] ));
        }
        mo.addprofile( new parameter(gdb_prof_str_K_DOMAIN_CENTERED, kc ) );
      }

    };
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - /
    static void checkprofiles( model_class & mo ){
      // we do a couple checks:
      // FCOEFFs need only 4 elements
      const vector<REAL> & fcoeffs = mo.get_global_prof_by_name(gdb_prof_str_FCOEFFS).getData();
      profcheck_assert(fcoeffs.size()>=4,"Profile '%s' needs at least '%i' elements.",gdb_prof_str_FCOEFFS,4);
      // K_DOMAIN must be monotonically increasing function with positive definite values
      const vector<REAL> & kdomain = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN).getData();
      profcheck_assert(kdomain.size()>=2,"Profile '%s' must have at least '%i' elements.",gdb_prof_str_K_DOMAIN,2);
      profcheck_assert(kdomain[0]>=0,"Profile '%s' must consist of positive definite elements.",gdb_prof_str_K_DOMAIN);
      for(unsigned int i = 1; i < kdomain.size(); i++){
        profcheck_assert(kdomain[i]>kdomain[i-1],"Profile '%s' must increase monotonically.",gdb_prof_str_K_DOMAIN);
      }
      if(!mo.does_prof_exist(gdb_prof_str_K_DOMAIN_CENTERED)){
        mapkcenteredprofile(mo);
      };
      if(!mo.does_prof_exist(gdb_prof_str_DK_DOMAIN)){
        mapdkprofile(mo); 
      }
      if( ! mo.does_prof_exist(gdb_prof_str_K_DOMAIN_SPACING) ){
        mapkspacingprofile(mo);
      }
      const vector<REAL> & dkdomain = mo.get_global_prof_by_name(gdb_prof_str_DK_DOMAIN).getData();
      profcheck_assert( dkdomain.size() == kdomain.size() -1 , 
          "Profile '%s' has size inconsistent with profile '%s'.", gdb_prof_str_DK_DOMAIN, gdb_prof_str_K_DOMAIN);
      // make sure all elements are greater than 0, warn if they do not match kdomain diffs
      for(size_t i = 0; i < dkdomain.size(); i++){
        const double dk_val = dkdomain[i];
        const double kn = kdomain[i+1];
        const double kn_1 = kdomain[i];
        profcheck_assert( dk_val > 0, "Profile '%s' must must consist of positive elements.", 
            gdb_prof_str_DK_DOMAIN);
        double diff = kn-kn_1;
        if( abs( (dk_val - diff) / dk_val ) > 0.0000000001 ){ // print a warning
          fprintf(stderr, "Warning: Element '%i' of profile '%s' do not agree with profile '%s'.\n",
              static_cast<int>(i+1), gdb_prof_str_DK_DOMAIN, gdb_prof_str_K_DOMAIN);
          fprintf(stderr, "         kn-kn-1 = %f but dk_val = %f.\n", diff, dk_val );
        }
      };
      // check k_domain_centered, must lie within each range of kdomain 
      // and must have number elements 1 less than kdomain
      const vector<REAL> & kdomain_centered = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_CENTERED).getData();
      profcheck_assert(
          kdomain_centered.size()==kdomain.size()-1,
          "Profile '%s' must have size %i less than profile '%s'.", 
          gdb_prof_str_K_DOMAIN_CENTERED, 1, gdb_prof_str_K_DOMAIN
          );
      for(size_t i = 1; i< kdomain.size(); i++){
        profcheck_assert( kdomain_centered[i-1] <= kdomain[i] && kdomain_centered[i-1] >= kdomain[i-1],
            "Elements of profile '%s' must lie between elements of profile '%s'.",
            gdb_prof_str_K_DOMAIN_CENTERED,
            gdb_prof_str_K_DOMAIN );
      }

      // add the K_DOMAIN_SQUARED stuff
      profcheck_assert( !mo.does_prof_exist(gdb_prof_str_K_DOMAIN_CENTERED_SQUARED), 
          "Profile '%s' is reserved and cannot be defined externally.", gdb_prof_str_K_DOMAIN_CENTERED_SQUARED);
      vector<REAL> k2(kdomain_centered.size(), 0);
      for(unsigned int i = 0; i < k2.size(); i++){
        k2[i] = kdomain_centered[i]*kdomain_centered[i];
      }
      mo.addprofile( new parameter(gdb_prof_str_K_DOMAIN_CENTERED_SQUARED, k2) );

      // check nkdomain profiles
      const vector<int> & nkdomain = mo.get_global_prof_by_name(gdb_prof_str_NK_DOMAIN).getiData();
      profcheck_assert(nkdomain.size() > 0 , "Profile '%s' must have at least 1 element.", gdb_prof_str_NK_DOMAIN);
      profcheck_assert(nkdomain[0] >= 0, "Profile '%s' must consist of positive definate elements.",
          gdb_prof_str_NK_DOMAIN);
      for(unsigned int i = 1; i < nkdomain.size(); i++){
        profcheck_assert(nkdomain[i] > nkdomain[i-1], "Profile '%s' must increase monotonically.",
            gdb_prof_str_NK_DOMAIN);
      }
      // check the kspacing profile
      const vector<REAL> kspacing = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_SPACING).getData();
      profcheck_assert(kspacing.size()==kdomain_centered.size(), "Profile '%s' must have size equal to profile '%s'.",
          gdb_prof_str_K_DOMAIN_SPACING, gdb_prof_str_K_DOMAIN_CENTERED);
      for(size_t i = 0; i < kspacing.size(); i++){
        profcheck_assert( kspacing[i] == gdb_enum_UNIFORMSPACING || kspacing[i] == gdb_enum_LOGSPACING , 
            "Profile '%s' takes on only two values: %s, %s.",
            gdb_prof_str_K_DOMAIN_SPACING, "1 --> uniform spacing, 2 --> logarithmic spacing.");
        if(kdomain[i] == 0.0){
          profcheck_assert( kspacing[i] != gdb_enum_LOGSPACING , "Profile '%s' cannot take on log value if lower bound on profile '%s' is zero.",
              gdb_prof_str_K_DOMAIN_SPACING, gdb_prof_str_K_DOMAIN);
        }
      }
      // make sure the size of all profiles is less than max(INDEX)  - 1 so that overflows do not occur
      for(size_t i = 0; i < mo.get_num_global_profs(); i++){
        const parameter & prof = mo.get_global_prof_by_index(i);
        size_t size = prof.GetSize();
        profcheck_assert( size < static_cast<size_t>(std::numeric_limits< INDEX >::max() ) ,
            "Profile '%s' has size that exceeds allowable INDEX size.", prof.GetName() );
      }
    };
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - /
    static void assigndefaultparams( model_class & mo ){
      mo.addparameter(new parameter( gdb_str_MODEL_NAME, static_cast<const char *>(gdb_param_MODEL_NAME)));
      mo.addparameter(new parameter( gdb_str_PI, static_cast<REAL>(gdb_param_PI)));
      mo.addparameter(new parameter( gdb_str_FOURIER_TRANSFORM_COEFFICIENT, static_cast<REAL>(gdb_param_FOURIER_TRANSFORM_COEFFICIENT)));
      mo.addparameter(new parameter( gdb_str_PQDOMAINOPTION, static_cast<const char *>(gdb_param_PQDOMAINOPTION)));
      mo.addparameter(new parameter( gdb_str_PQFINEDOMAINACCURACYTHRESHOLD, static_cast<REAL>(gdb_param_PQFINEDOMAINACCURACYTHRESHOLD)));
      mo.addparameter(new parameter( gdb_str_PQFINEDOMAINSTARTINGRESOLUTION, static_cast<int>(gdb_param_PQFINEDOMAINSTARTINGRESOLUTION)));
      mo.addparameter(new parameter( gdb_str_PQFINEDOMAINRESOLUTIONFACTOR, static_cast<int>(gdb_param_PQFINEDOMAINRESOLUTIONFACTOR)));
      mo.addparameter(new parameter( gdb_str_NUMTHREADS, static_cast<int>(-1)));
      mo.addparameter(new parameter( gdb_str_NUMFIELDS, static_cast<int>(gdb_param_NUMFIELDS)));
      mo.addparameter(new parameter( gdb_str_DEBUGONERROR, static_cast<int>(gdb_param_DEBUGONERROR)));
      mo.addparameter(new parameter( gdb_str_TIMESTEP, static_cast<REAL>(gdb_param_TIMESTEP)));
      mo.addparameter(new parameter( gdb_str_IMPLICIT_DIAG_VARS, static_cast<int>(gdb_param_IMPLICIT_DIAG_VARS)));
      mo.addparameter(new parameter( gdb_str_INITIAL_NESTING_DEPTH, static_cast<int>(gdb_param_INITIAL_NESTING_DEPTH)));
      mo.addparameter(new parameter( gdb_str_MODEL_MESSAGE_OUTPUT_FILE, static_cast<const char *>(gdb_param_MODEL_MESSAGE_OUTPUT_FILE)));
      mo.addparameter(new parameter( gdb_str_MODEL_OUTPUT_FILE_TEMPLATE, static_cast<const char *>(gdb_param_MODEL_OUTPUT_FILE_TEMPLATE)));
      mo.addparameter(new parameter( gdb_str_MODEL_FIFO, static_cast<const char *>(gdb_param_MODEL_FIFO) ) );
      mo.addparameter(new parameter( gdb_str_NONLINEAR, static_cast<int>(gdb_param_NONLINEAR)));
      mo.addparameter(new parameter( gdb_str_STATIC_COEFFS_FILE, gdb_param_STATIC_COEFFS_FILE) );
      mo.addparameter(new parameter( gdb_str_INCLUDE_SEARCH_DIR, static_cast<const char *>(gdb_param_INCLUDE_SEARCH_DIR) ) );
      mo.addparameter(new parameter( gdb_str_INCLUDE_SEARCH_DIR_ENV_VAR_NAME, static_cast<const char *>(gdb_param_INCLUDE_SEARCH_DIR_ENV_VAR_NAME) ) );
      mo.addparameter(new parameter( gdb_str_IGNORESLOT, static_cast<int>(gdb_param_IGNORESLOT) ));
      mo.addparameter(new parameter( gdb_str_SLOT_WIDTH_SCALE, static_cast<REAL>(gdb_param_SLOT_WIDTH_SCALE)));
      mo.addparameter(new parameter( gdb_str_IGNORETRIANGLEBASE, static_cast<int>(gdb_param_IGNORETRIANGLEBASE) ));
      mo.addparameter(new parameter( gdb_str_KCUTOFFMULTIPLIER, static_cast<REAL>(gdb_param_KCUTOFFMULTIPLIER) ) );
      mo.addparameter(new parameter( gdb_str_IGNOREREALIZABILITYCONSTRAINTS, static_cast<int>(gdb_param_IGNOREREALIZABILITYCONSTRAINTS) ) );

      mo.addparameter(new parameter( gdb_str_ATMOSPHERE_LAYER_1_FIELD_NUMBER, static_cast<int>(gdb_param_ATMOSPHERE_LAYER_1_FIELD_NUMBER) ) );
      mo.addparameter(new parameter( gdb_str_ATMOSPHERE_LAYER_2_FIELD_NUMBER, static_cast<int>(gdb_param_ATMOSPHERE_LAYER_2_FIELD_NUMBER) ) );
      mo.addparameter(new parameter( gdb_str_OCEAN_LAYER_1_FIELD_NUMBER, static_cast<int>(gdb_param_OCEAN_LAYER_1_FIELD_NUMBER) ) );
      mo.addparameter(new parameter( gdb_str_OCEAN_LAYER_2_FIELD_NUMBER, static_cast<int>(gdb_param_OCEAN_LAYER_2_FIELD_NUMBER) ) );
      mo.addparameter(new parameter( gdb_str_MIXED_LAYER_FIELD_NUMBER, static_cast<int>(gdb_param_MIXED_LAYER_FIELD_NUMBER) ) );
      mo.addparameter(new parameter( gdb_str_TOPO_LAYER_FIELD_NUMBER, static_cast<int>(gdb_param_TOPO_LAYER_FIELD_NUMBER) ) );
    };
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -/
    static void assigndefaultprofiles( model_class & mo ){
      //parameter * p = NULL;

      vector<REAL> kdomain;
      kdomain.push_back(1);
      kdomain.push_back(2);
      kdomain.push_back(3);
      kdomain.push_back(4);

      mo.addprofile(new parameter ( gdb_prof_str_K_DOMAIN, kdomain ) );

      vector<int> nkdomain;
      nkdomain.push_back(0);
      nkdomain.push_back(1);

      mo.addprofile( new parameter ( gdb_prof_str_NK_DOMAIN, nkdomain ) );

      vector<REAL> fcoeffs(4, 0.0 );

      mo.addprofile( new parameter ( gdb_prof_str_FCOEFFS, fcoeffs ) );

      vector<int> aasampfreq(1, 2);
      mo.addprofile( new parameter ( gdb_prof_str_AASAMPLINGFREQUENCY, aasampfreq ) );

      vector<REAL> multipliers;
      multipliers.push_back( 10.0 );
      multipliers.push_back( 50.0 );
      multipliers.push_back( 100.0 );
      multipliers.push_back( 500.0 );
      multipliers.push_back( 1000.0 );
      multipliers.push_back( 2000.0 );
      multipliers.push_back( 3000.0 );
      multipliers.push_back( 5000.0 );

      mo.addprofile( new parameter( gdb_prof_str_FINEKDOMAINMULTIPLIERS, multipliers ) );


    };
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -/
    static bool istokendefinedinparser(const ParserX & p, const string & token)
    {
      const var_maptype & vars = p.GetVar();

      var_maptype::const_iterator item = vars.begin();
      for(; item != vars.end(); ++item){
        if(item->first == token){ return true; }
      }

      // check functions
      const fun_multimap & fmap = p.GetFunDef();
      for( fun_multimap::const_iterator item = fmap.begin(); item != fmap.end(); ++item){
        if(item -> first == token ) return true;
      }

      // now, do the constants
      const val_maptype & cmap = p.GetConst();
      if (cmap.size())
      {
        val_maptype::const_iterator item = cmap.begin();
        for(; item != cmap.end(); ++item){
          if(item->first == token){ return true; }
        }
      }
      return false;
    }
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -//
    static bool save_model_profiles(const model_class & mo, NcFile * ncfile){
      const size_t num_profs = mo.get_num_global_profs();
      for(size_t i = 0; i < num_profs; i++){
        const parameter & prof = mo.get_global_prof_by_index(i);
        const parameter::type ptype = prof.getType();
        const char * name = prof.GetName();
        switch(ptype){
          case parameter::Real :
            { const vector<REAL> & data = prof.getData();
              REAL * array = new REAL[data.size()];
              for(size_t x = 0; x < data.size(); x++) array[x] = data[x];
              ncfile->add_att(name, data.size(), array);
              delete [] array;
            }
            break;
          case parameter::Integer :
            { const vector<int> & data = prof.getiData();
              int * array = new int[data.size()];
              for(size_t x = 0; x < data.size(); x++) array[x] = data[x];
              ncfile->add_att(name, data.size(), array);
              delete [] array;
            }
            break;
          case parameter::String :
            // each is added individually
            break;
        }; 
      }
      return true;
    };
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -//
    static bool save_model_parameters(const model_class & mo, NcFile * ncfile){
      NcError error_status( NcError::silent_nonfatal );
      const size_t num_params = mo.get_num_global_params();
      for(size_t i = 0; i < num_params; i++){
        const parameter & p = mo.get_global_param_by_index(i);
        const char * name = p.GetName();
        const parameter::type ptype = p.getType();
        switch(ptype){
          case parameter::Real :
            ncfile->add_att(name, p.getfValue()); 
            break;
          case parameter::String :
            ncfile->add_att(name, p.getsValue());
            break;
          case parameter::Integer :
            ncfile->add_att(name, p.getiValue());
            break;
        };
      }
      if( error_status.get_err() == NC_NOERR ) return true;
      mo.printfmsg("NETCDF Error: %s.\n", error_status.get_errmsg());
      return false;
    };
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -//
    static bool get_att_string(NcAtt * att, string & dest)
    {
      NcError error_status( NcError::silent_nonfatal );
      char * asstring = att->as_string(0);
      if(asstring){
        dest.assign(asstring);
        delete asstring;
        return true;
      }
      fprintf(stderr, "NETCDF Error: %s.\n", error_status.get_errmsg());
      fprintf(stderr, "Unable to read attribute from file.\n");
      return false;
    }
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -// 
#define declare_string_attribute(dest, att) string dest; { if(!get_att_string(att,dest)) return false; };
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -//
    static bool create_script_path_file_name(const char * stemplate, string & dest)
    { // this function is rather crude at this point, it just needs to generate
      // a file name that is determined by the incoming template (say out.nc) and 
      // the current time/date etc
      dest.assign(stemplate);
      dest+=".script.info.nc";
      return true;
    }
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -//
    static void get_function_lines(const preprocessordata & pdata, const int nid, const int sline, const int eline,string& dest){
      const namelistdata & ndata = *(pdata.namelists[nid]);  /// .getnamelistbyid(nid);
      const string & ncontents = ndata.filecontent;
      int ln = 1;
      int spos = -1;
      int epos = -1;
      for(size_t pos = 0; pos < ncontents.size(); ++pos ){
        if(ln >= sline && spos < 0) spos = pos;
        if(ln <= eline ) epos = pos;
        if(ln > eline) break;
        if( ncontents[pos] == '\n' ) ln++;
      }
      if(epos>=spos){
        dest = ncontents.substr(spos, epos-spos+1);
      }
    }
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -//
    static string getidcode(const char * stemplate)
    {
      // first generate number
      srand ( time(NULL) ); // init the seed
      int rNumb = rand() % 100000 + 1;
      string rvalue = string(stemplate) + "::" + convertint2string(rNumb);
      time_t rtime; 
      struct tm * timeinfo;
      time( & rtime );
      timeinfo = localtime( &rtime ); // returns statically allocate timeinfo object
      string date( asctime(  timeinfo ) ); // returns statically allocated string
      for(int i = 0; i < static_cast<int>(date.size()); i++){
        if(date[i] == ' ') date[i]='_'; // get rid of any spaces
        if(date[i] == '\n') date[i]=0; // get rid of any spaces
        if(date[i] == ':') date[i]='-'; // get rid of any spaces
      }
      rvalue += "::" + date;
      return rvalue;
    }
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -//
    static bool save_model_script_info(const model_class & mo, NcFile * ncfile, const char * template_name){
      NcError error_status( NcError::silent_nonfatal );

      vector< const model_diagnostics_block * > func_list;

      mo.getfunctiondiaglist( func_list ); // retrieve pointers to diag blocks that are functions

      string script_path;
      if(!create_script_path_file_name(template_name, script_path)){
        mo.printfmsg("Error: Could not open %s for writing script information.\n");
        return false;
      }
      remove(script_path.c_str());

      /* write the output data to the netcdf file */
      NcVar * var = ncfile->add_var("ScriptInformation", ncInt, NULL);
      if( var == NULL || error_status.get_err() != NC_NOERR )
      {
        mo.printfmsg("NETCDF Error: %s.\n", error_status.get_errmsg());
        mo.printfmsg("Failed writing data pointing to location of script information.\n");
        return false;
      }
      var -> add_att("path", script_path.c_str() );
      string idcode = getidcode(template_name);
      var -> add_att("idcode", idcode.c_str());

      /* attempt to open the new file as a netcdf file */
      nchandler nc_container;
      NcFile::FileMode fmode = NcFile::New;
      double nchandle = 0;
      nchandle = nc_container.createncfile(script_path.c_str(), fmode);
      NcFile * scriptfile = nc_container.getncfile(static_cast<int>(nchandle));
      scriptfile -> add_att("idcode", idcode.c_str());

      if(nchandle == 0){ // error occured
        mo.printfmsg("File '%s' could not be created for writing sript information.\n", script_path.c_str());
        return false;
      }

      for(size_t iter = 0; iter < func_list.size(); iter++){ // iterate through all the functions, etc, write data 
        const model_diagnostics_block & func = *(func_list[iter]);
        const string & name = func.getfuncname();
        string nc_name("script_block_"); nc_name+=convertint2string(iter);
        NcVar * var = scriptfile->add_var(nc_name.c_str(), ncChar, NULL);
        if(var == NULL || error_status.get_err() != NC_NOERR){
          mo.printfmsg("NETCDF Error: %s.\n", error_status.get_errmsg());
          mo.printfmsg("Encountered error while writing script info.\n");
          return false;
        }
        var->add_att("name", name.c_str() );
        if( func.isblockfunction()){
          const int numargs = func.getfuncnumargs();
          const vector<char_type> cstacktypes = func.GetCallStackTypes(); 
          const char_type rtype = func.GetReturnType();
          string args;
          create_function_call_stack_string( cstacktypes, args );
          var->add_att("type", "function");
          var->add_att("callsignature", args.c_str());
          var->add_att("arguments", numargs);
          var->add_att("returns", get_type_id(rtype) );
        }else{
          var->add_att("type", "script");
        }
        string data;
        //const vector<model_diagnostics_block::instruction *> & instructions = func.instructions;
        int nid = func.debug_attributes.nid;
        int sline =func.debug_attributes.slinenumber;
        int eline = func.debug_attributes.elinenumber;
        // grab the string 
        get_function_lines( mo.preprocdata, nid, sline, eline, data );
        var->add_att("body", data.c_str() );
      }

      if( error_status.get_err() == NC_NOERR ) return true;
      mo.printfmsg("NETCDF Error: %s.\n", error_status.get_errmsg());
      return false;
    };
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -//
    static const char * nctype2string(NcType nctype){
      switch(nctype){
        case ncInt:
          return "INTEGER";
        case ncDouble:
        case ncFloat:
          return "REAL";
        case ncChar:
          return "STRING";
        default :
          return "unknown";
      };
    };
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -//
    static const char * ptype2string(parameter::type ptype){
      switch(ptype){
        case parameter::Real:
          return "REAL";
        case parameter::Integer:
          return "INTEGER";
        case parameter::String:
          return "STRING";
        default:
          return "unknown";
      }
    };
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -//
    static bool restore_check_param_profiles(const model_class & mo, NcFile * ncfile){
      // all parameters and profiles loaded in global attributes must agree
      // in value with those in the model !
      const int num_atts = ncfile->num_atts();
      for(int iter = 0; iter < num_atts; iter++){
        NcAtt * att = ncfile->get_att(iter); // get the nth attribute
        const char * a_name = att->name();
        if( mo.does_param_exist(a_name) ){
          // this is a parameter, compare values
          const parameter & param = mo.get_global_param_by_name( a_name );
          NcType vtype = att->type();
          // make sure sizes are 1 for attribute
          const int asize = static_cast<int>(att->num_vals());
          if( vtype != ncChar && asize > 1 ){
            fprintf(stderr, "Parameter '%s' mismatch with current model instantiation.\n", a_name);
            fprintf(stderr, "Saved state has more than 1 element (i.e. profile) but model expects a parameter.\n"); 
            return false;
          }
          parameter::type ptype = param.getType();
          if(vtype == ncInt && ptype == parameter::Integer){
            int p_value = param.getiValue();
            int a_value = att->as_int(0);
            if( a_value != p_value ){
              fprintf(stderr, "Warning: Parameter '%s' mismatch with current model instantiation.\n", a_name);
              fprintf(stderr, "Warning: Model has defined value of '%i' but saved state has a value of '%i'.\n", p_value, a_value);
            }
          }else if( ptype == parameter::Real && ( vtype == ncFloat || vtype == ncDouble || vtype == ncInt ) ){
            float p_value = static_cast<float>( param.getfValue() );
            float a_value = att->as_float(0);
            if( a_value != p_value ){
              fprintf(stderr, "Warning: Parameter '%s' mismatch with current model instantiation.\n", a_name);
              fprintf(stderr, "Warning: Model has defined value of '%f' but saved state has a value of '%f'.\n", p_value, a_value);
            }
          }else if(vtype == ncChar && ptype == parameter::String){
            const string p_value( param.getsValue() );
            declare_string_attribute(a_value,  att);
            if ( a_value != p_value ){
              fprintf(stderr, "Warning: Parameter '%s' mismatch with current model instantiation.\n", a_name);
              fprintf(stderr, "Warning: Model has defined value of \"%s\" but saved state has a value of \"%s\".\n", 
                  p_value.c_str(), a_value.c_str());
            }
          }else{
            // error type mismatch
            fprintf(stderr, "Parameter '%s' type mismatch with current model instantiation.\n", a_name);
            fprintf(stderr, "Model's version of the parameter has type '%s' but saved state has type '%s'.\n", nctype2string(vtype),
                ptype2string(ptype));
            return false;
          }
        }else if( mo.does_prof_exist( a_name ) ){
          const parameter & prof = mo.get_global_prof_by_name( a_name );
          NcType vtype = att->type();
          parameter::type ptype = prof.getType();
          // make sure size is the same
          const int psize = static_cast<int>(prof.GetSize());
          const int asize = static_cast<int>(att->num_vals());
          if( psize != asize ){
            fprintf(stderr, "Profile '%s' size mismatch with current model instantiation.\n", a_name);
            fprintf(stderr, "Model's version has size '%i' but version in the saved state has size '%i'.\n", psize, asize);
            return false;
          }
          for(int iter = 0; iter < psize; iter++){
            if(vtype == ncInt && ptype == parameter::Integer){
              int p_value = prof.getiData()[iter];
              int a_value = att->as_int(iter);
              if( a_value != p_value ){
                fprintf(stderr, "Warning: Profile '%s' mismatch with current model instantiation in element '%i'.\n", a_name, iter);
                fprintf(stderr, "Warning: Model has defined value of '%i' but saved state has a value of '%i'.\n", p_value, a_value);
              }
            }else if( ptype == parameter::Real && ( vtype == ncFloat || vtype == ncDouble || vtype == ncInt ) ){
              float p_value = static_cast<float>(prof.getData()[iter]);
              float a_value = att->as_float(iter);
              if( a_value != p_value ){
                fprintf(stderr, "Warning: Profile '%s' mismatch with current model instantiation in element '%i'.\n", a_name, iter);
                fprintf(stderr, "Warning: Model has defined value of '%f' but saved state has a value of '%f'.\n", p_value, a_value);
                return false;
              };
            }else if( ptype == parameter::String ){ // just ignore this prof

            }else{ // bad type or mismatch 
              fprintf(stderr, "Warning: Profile '%s' type mismatch with current model instantiation.\n", a_name);
              fprintf(stderr, "Warning: Model's version of the profile has type '%s' but saved state has type '%s'.\n", nctype2string(vtype),
                  ptype2string(ptype));
              return false;
            }
          }
        }else{ // uh oh parameter/prof is missing
          if( att->num_vals() > 1 ){
            fprintf(stderr, "Found parameter '%s' in saved state, but parameter does not exist in current instantiation of the model.\n",
                a_name);
          }else{
            fprintf(stderr, "Found profile '%s' in saved state, but parameter does not exist in current instantiation of the model.\n",
                a_name);
          }
          return false;
        }
      }
      return true;
    };
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -//
    static bool print_original_function_body(NcVar * func){
      NcError error_status( NcError::silent_nonfatal );
      declare_string_attribute(fname, (func->get_att("name")));
      fprintf(stderr, " - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
      fprintf(stderr, "Function '%s' as originally stored in the model's saved state:\n", fname.c_str());
      fprintf(stderr, " - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
      NcAtt * obody_att = func->get_att("body");
      declare_string_attribute(body, obody_att);
      fprintf(stderr, "%s", body.c_str());
      fprintf(stderr, " - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
      fprintf(stderr, " Note: Only the number of arguments and the function name must agree. \n");
      fprintf(stderr, "       The body of the function, its argument and return types can differ\n");
      fprintf(stderr, "       between the current instantiation of the model and the saved state.\n");
      fprintf(stderr, " - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
      return true;
    };
    //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -//
    static bool restore_check_script_functions(const model_class & mo, NcFile * ncfile){
      NcError error_status( NcError::silent_nonfatal );

      nchandler nc_container;
      vector< const model_diagnostics_block * > func_list;
      mo.getfunctiondiaglist( func_list ); // retrieve pointers to diag blocks that are functions

      // grab file path
      NcVar * script_var = ncfile->get_var("ScriptInformation");
      if(!script_var){ fprintf(stderr, "Error: Could not extract script information.\n"); return false; }
      NcAtt * spath_att = script_var->get_att("path");
      NcAtt * idcode_att= script_var->get_att("idcode");
      if(!spath_att || !idcode_att){ fprintf(stderr, "Error: Could not extract script information.\n"); return false; }
      declare_string_attribute(spath, spath_att);
      declare_string_attribute(idcode, idcode_att);

      // now open the script information file
      double nchandle = 0;
      // the file must exist, check to make sure it does 
      nchandle = nc_container.createncfile(spath.c_str(), NcFile::ReadOnly);
      if(nchandle == 0){
        fprintf(stderr, "Could not open auxiliary file '%s' that contains the model state's script information.\n", spath.c_str());
        return false;
      }
      NcFile * script_file = nc_container.getncfile(static_cast<int>(nchandle));
      // check idcode
      NcAtt * spt_idcode_att = script_file->get_att("idcode");
      if( ! spt_idcode_att ){ fprintf(stderr, "Could not retrieve idcode for file '%s'.\n", spath.c_str()); return false; }
      declare_string_attribute( spt_idcode, spt_idcode_att );
      if( idcode != spt_idcode ){
        fprintf(stderr, "Mismatched idcode in and auxiliary file '%s'. Corrupted data?\n", spath.c_str()); return false; 
      }

      int iter = 0;
      while(true){
        string nc_name("script_block_"); nc_name += convertint2string(iter);
        NcVar * var = script_file->get_var(nc_name.c_str());
        if( var == NULL ) break; // otherwise peruse the block
        NcAtt * fname_att = var->get_att("name");
        NcAtt * ntype_att = var->get_att("type");  
        if(fname_att == NULL || ntype_att == NULL ){
          fprintf(stderr, "Error: Could not extract information for script function #%i from saved state.\n", iter);
          fprintf(stderr, "       File may be corrupt.\n");
          return false;
        }
        declare_string_attribute(afname, fname_att);
        declare_string_attribute(stype, ntype_att);
        if(stype != "function"){
          iter++; continue;
        }
        NcAtt * fsize_att = var->get_att("arguments");
        const int afsize = fsize_att->as_int(0);

        // now look up function in current model instantiation and make sure call stack agrees
        bool found_function = false;
        for(size_t u = 0; u < func_list.size(); u++){
          if( func_list[u]->getfuncname() == afname){
            bool mismatch = (afsize != func_list[u]->getfuncnumargs());
            if(mismatch){ // number or arguments must agree, types not necessarily so, print a warning if types differ
              fprintf(stderr, "Error: Mismatch in function '%s' definition with current instantiation of the model.\n", afname.c_str());
              fprintf(stderr, "The model has a function definition that takes '%i' arguments, but the saved state takes '%i' arguments.\n",
                  func_list[u]->getfuncnumargs(), afsize);
              fprintf(stderr, "Function call stacks must agree otherwise compatability between current model\n");
              fprintf(stderr, "instantiation and saved state cannot be gauranteed.\n");
              print_original_function_body(var); 
              return false;
            }  
            declare_string_attribute(rtypestr, var->get_att("returns"));
            mismatch = ( rtypestr != get_type_id( func_list[u]->GetReturnType() ) );
            // now check the call stack
            const vector<char_type> cstacktypes = func_list[u]->GetCallStackTypes(); 
            string args;
            create_function_call_stack_string( cstacktypes, args );    // check entire call stack next
            declare_string_attribute( saved_args, var->get_att("callsignature"));
            mismatch = mismatch || ( saved_args != args );
            if( mismatch ){
              fprintf(stderr, "Warning: Mismatch in function '%s' definition with current instantiation of the model.\n", afname.c_str());
              fprintf(stderr, "Current instantiation of the model has a function defined as '%s %s%s'\n",
                  get_type_id(func_list[u]->GetReturnType()), afname.c_str(), args.c_str());
              fprintf(stderr,"\tbut saved state is defined as '%s %s%s'\n",rtypestr.c_str(), afname.c_str(), saved_args.c_str());
              fprintf(stderr, "Function call stacks should agree on types otherwise compatability between current model\n");
              fprintf(stderr, "instantiation and saved state cannot be guaranteed.\n");
              print_original_function_body(var); 
            } 
            found_function = true;
            break;
          }
        }
        if( ! found_function ){ // error
          fprintf(stderr, "Error: Function '%s' is missing from current model instantiation.\n", afname.c_str());
          fprintf(stderr, "The saved state requires the definition of this function with '%i' arguments.\n", afsize);
          fprintf(stderr, "All user defined functions present in the saved state must exist in current model instantiation\n");
          fprintf(stderr, "otherwise compatability between current model and the saved state\n");
          fprintf(stderr, "cannot be guaranteed.\n");
          print_original_function_body(var); 
          return false;
        }   
        ++iter;
      }
      return true;
    };
}; // end model_class::local
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
bool istokendefinedinparser(const ParserX & p, const string & token)
{
  return model_class::local::istokendefinedinparser(p, token);
}
/**************************************************************************************
  class model_class
 **************************************************************************************/
model_class::model_class() : preprocdata(this), mdebugger(this)
{
  this->model_state.preprocessed = false;
  this->model_state.compiled = false;
  this->model_state.io_initialized = false;
  this->model_state.physics_initialized = false;
  this->model_state.state_initialized = false;

  this->default_labels.namelist.assign("default.qg");
  this->default_labels.stdout_name.assign("");
  this->default_labels.stdin_name.assign("");
  this->default_labels.stdmsg_name.assign("");

  // set default field values
  this->atmosphere_layer_1 = 1;
  this->atmosphere_layer_2 = 2;
  this->ocean_layer_1 = 3;
  this->ocean_layer_2 = 4;
  this->mixed_layer = 5;
  this->topo_layer = 6;
};
//-------------------------------------------------------------------//
bool model_class::init_input(const string & label){
  const char * ifname = label.size() > 0 ? label.c_str() : NULL;
  // default value is to set stdin to stdin
  if( !ifname ){
    if(this->does_param_exist(gdb_str_STDIN)){
      if( this->get_global_param_by_name(gdb_str_STDIN).getType() == parameter::String){
        ifname = this->get_global_param_by_name(gdb_str_STDIN).getsValue();
      }
    }
  } 
  if(ifname ){
    if( ! this->STDIN.open(ifname, "r") ){
      // FILE * input = fopen(ifname, "a");
      fprintf(stderr, "Cannot open file '%s' for reading.\n", ifname);
      fprintf(stderr, "Unable to initialize stdin.\n");
      return false;
    }
    assert( ! (this->STDIN.is_closed() ) );
  }else{
    this->STDIN.open(stdin); // default behavior
  }
  return true;
};
//-------------------------------------------------------------------//
int model_class::get_num_threads() const{
#ifdef SERIALQGCM
  return 1;
#else
  int num_procs_var = this->get_global_param_by_name(gdb_str_NUMTHREADS).getiValue();
  if(num_procs_var > 0){
    return num_procs_var;
  }else{
    // return the number of processors available on the machine
    return omp_get_num_procs(); 
  }
#endif
};
//-------------------------------------------------------------------//
bool model_class::init_master_output(const string & lstd, const string & lmsg){
  const char * stdname = lstd.size() > 0 ? lstd.c_str() : NULL;
  if(!stdname){ // try loading data from a parameter
    if(this->does_param_exist(gdb_str_STDOUT)){
      const parameter & pstdout = this->get_global_param_by_name(gdb_str_STDOUT);
      if(pstdout.getType() == parameter::String){
        stdname = pstdout.getsValue();
      }
    }
  }
  // allocate handles
  if(!stdname){
    (this->STDOUTs[0]).open(stdout);
    //this->STDOUTs[0] = stdout; // default to stdout
  }else{
    // Make sure the file does not already exist
    const char * errormsg = NULL;
    if( ! check_stream_output_status(stdname, errormsg) ){
      fprintf(stderr, "Error opening '%s' for writing model output.\n", stdname);
      fprintf(stderr, "Reason: %s\n", errormsg);
      return false;
    }
    if( ! this->STDOUTs[0].open( stdname, "w" ) ){
      //this->STDOUTs[0] = fopen(stdname, "w");
      //if( !this->STDOUTs[0] ){
      fprintf(stderr, "Could not open file '%s' for writing model output.\n", stdname);
      return false;
    };
    }
    const char * msgname = lmsg.size() > 0 ? lmsg.c_str() : NULL;
    // now handle the messages part
#ifndef SERIALQGCM 
    if(!msgname){ // use the default model.msg file for message output 
      msgname = this->get_global_param_by_name(gdb_str_MODEL_MESSAGE_OUTPUT_FILE).getsValue();
    }
#endif
    const char * errormsg = NULL;
    if( ! check_stream_output_status(msgname, errormsg) ){
      fprintf(stderr, "Error opening '%s' for writing model messages.\n", msgname);
      fprintf(stderr, "Reason: %s\n", errormsg);
      return false;
    };
    if(msgname){
      //this->MSGOUTs[0] = fopen(msgname, "w");
      if(! this->MSGOUTs[0].open(msgname, "w")){
        //if( ! this->MSGOUTs[0] ){
        fprintf(stderr, "Could not open file '%s' for writing model messages.\n", msgname);
        return false;
      } 
      }else{
        // this->MSGOUTs[0] = stderr;
        this->MSGOUTs[0].open(stderr);
      };
      return true; 
    };
#ifdef SERIALQGCM
    //-------------------------------------------------------------------//
    bool model_class::init_output(const string & stdname, const string & msgname){
      this->out_handle_size = 1;
      // this->STDOUTs = new FILE * [this->out_handle_size];
      // this->MSGOUTs = new FILE * [this->out_handle_size];
      this->STDOUTs.resize( this->out_handle_size, qstream() );
      this->MSGOUTs.resize( this->out_handle_size, qstream() );
      if( ! init_master_output(stdname, msgname) ) return false;
      if( this->getSTDMSG().getpath().size() > 0 ){
        if( !this->mdebugger.setDBGOUT( this->getSTDMSG().getpath().c_str() )) return false;
        if( !this->mdebugger.setDBGIN( this->getSTDIN().getpath().c_str() )) return false;
      }else{
        if( ! this->mdebugger.setDBGOUT(stderr) ) return false;
        if( ! this->mdebugger.setDBGIN(stdin) ) return false;
      }
      return true; 
    };
#else
    //-------------------------------------------------------------------//
    bool model_class::init_output(const string & stdname, const string & msgname){
      const size_t NUMTHREADS = this->get_num_threads();
      this->out_handle_size = NUMTHREADS + 1;
      rassert( this->out_handle_size >= 1 );
      this->STDOUTs.resize( this->out_handle_size, qstream() );
      this->MSGOUTs.resize( this->out_handle_size, qstream() );
      if( ! init_master_output(stdname, msgname) ) return false;
      const char * msgtemplate = this->get_global_param_by_name(gdb_str_MODEL_MESSAGE_OUTPUT_FILE).getsValue();
      const char * stdtemplate = this->get_global_param_by_name(gdb_str_MODEL_OUTPUT_FILE_TEMPLATE).getsValue();
      for(int i = 1; i < this->out_handle_size; i++){
        string strmsgname(msgtemplate);
        string strstdname(stdtemplate);
        string number = string(1,i+'0');
        strmsgname = strmsgname + string(".") + number;
        strstdname = strstdname + string(".") + number;
        if(! this->MSGOUTs[i].open(strmsgname.c_str(), "w")){
          fprintf(stderr, "File '%s' cannot be opened for writing model messages.\n", strmsgname.c_str());
          return false;
        }
        if( ! this->STDOUTs[i].open( strstdname.c_str(), "w" ) ){
          fprintf(stderr, "File '%s' cannot be opened for writing model stdout.\n", strstdname.c_str());
          return false;
        }
      }
      return true;
    };
#endif
    //-------------------------------------------------------------------//
    void model_class::preprocess(){
      if( ! this->model_state.preprocessed ){
        const char * namelist = this->default_labels.namelist.c_str();
        local::assigndefaultparams( *this );
        local::assigndefaultprofiles( * this );
        preprocessor::invoke(this->preprocdata, namelist);
        this->model_state.preprocessed = true;
      }
    };
    //-------------------------------------------------------------------//
    int model_class::runscript(const vector<const char *> & scriptlist){
      if( this->model_state.compiled ) return -1;
      this->preprocess();
      // now load the requested file
      for(size_t u = 0; u < scriptlist.size(); u++){
        this->preprocdata.readnamelist(scriptlist[u]);
      }
      // now compile the model
      this->compile();
      this->init_io();
      for(size_t i = 0; i < mdiags.size(); i++){
        if( mdiags[i]->isblockscript() ){
          // set args
          ptr_val_type * args = new ptr_val_type[1];
          Value arg;
          args[0].Reset(&arg);
          // args[0] = 0; // may support an array argument later
          mdiags[i]->assignarglist(args, 1 );
          bool rvalue = mdiags[i]->runscriptblock(debugonerror());
          // delete [] args; (deconstructors called automatically)
          if( rvalue != 0 ) return static_cast<int>(rvalue);
        }
      }
      return 0;
    };

    //-------------------------------------------------------------------//
    void model_class::init_io(){ //const char * in_name, const char * stdout_name, const char * msg_name){
      if( ! this->model_state.compiled ) this->compile();
      if( ! this->model_state.io_initialized ){
        const char * in_name = this->default_labels.stdin_name.c_str();
        const char * stdout_name = this->default_labels.stdout_name.c_str();
        const char * msg_name = this->default_labels.stdmsg_name.c_str();
        if(!this->init_input(in_name)){
          throw cgcmexception("Could not initialize model standard input.", -1);
        }
        if(!this->init_output(stdout_name, msg_name)){ // initialize serial and/or parallelized output
          throw cgcmexception("Could not initialize model output.", -1);
        };
        this->model_state.io_initialized = true;
      }
    };
    //-------------------------------------------------------------------//
    void model_class::init_state(){ // create S and P system and initialize to model for integration
      if( ! this->model_state.compiled ) this->compile();
      if( ! this->model_state.io_initialized ) this->init_io();
      if( ! this->model_state.field_geometry_initialized ) this->init_field_maps();
      if( ! this->model_state.physics_initialized ) this->init_physics();
      if( ! this->model_state.state_initialized ){
        this->smodule.init(*this, this->spre_data);
        this->pmodule.init(this);
        this->model_state.state_initialized = true;
      };
    };
    //-------------------------------------------------------------------//
    void model_class::init_field_maps(){
      if( ! this->model_state.compiled ) this->compile();
      if( ! this->model_state.io_initialized ) this->init_io();
      if( ! this->model_state.field_geometry_initialized ){
        this->atmosphere_layer_1 = static_cast<INDEX>(this->get_global_param_by_name(gdb_str_ATMOSPHERE_LAYER_1_FIELD_NUMBER).getiValue());
        this->atmosphere_layer_2 = static_cast<INDEX>(this->get_global_param_by_name(gdb_str_ATMOSPHERE_LAYER_2_FIELD_NUMBER).getiValue());
        this->ocean_layer_1 = static_cast<INDEX>(this->get_global_param_by_name(gdb_str_OCEAN_LAYER_1_FIELD_NUMBER).getiValue());
        this->ocean_layer_2 = static_cast<INDEX>(this->get_global_param_by_name(gdb_str_OCEAN_LAYER_2_FIELD_NUMBER).getiValue());
        this->mixed_layer = static_cast<INDEX>(this->get_global_param_by_name(gdb_str_MIXED_LAYER_FIELD_NUMBER).getiValue());
        this->topo_layer = static_cast<INDEX>(this->get_global_param_by_name(gdb_str_TOPO_LAYER_FIELD_NUMBER).getiValue());
        this->model_state.field_geometry_initialized = true;
      }
    };
    //-------------------------------------------------------------------//
    void model_class::init_physics(){
      if( ! this->model_state.compiled ) this->compile();
      if( ! this->model_state.physics_initialized ){
        // assign the field numbers if needed
        this->init_field_maps();

        bool nonlinear = (*this).get_global_param_by_name(gdb_str_NONLINEAR).getiValue() != 0;
        const char * fname = (*this).get_global_param_by_name(gdb_str_STATIC_COEFFS_FILE).getsValue();
        nchandler nc_container;
        double nchandle;
        if( nonlinear ){
          if(! mx_read_static_data( *this, fname, this->rcoeffs, this->scoeffs, this->Acoeffs, this->Qcoeffs,this->QIcoeffs))
          {
            throw cgcmexception(string("Failed to read static coefficients file: '") + string(fname) + string("'."), -15);
          }
        }else{
          if( ! mx_read_static_data( *this, fname, this->Qcoeffs, this->QIcoeffs ) ){
            throw cgcmexception(string("Failed to read static coefficients file: '") + string(fname) + string("'."), -15);
          }
        }
        this->model_state.physics_initialized = true;
      }
    };
    //-------------------------------------------------------------------//
    const string model_class::get_default_stdin() const{
      return this->default_labels.stdin_name.c_str();
    };
    //-------------------------------------------------------------------//
    const string model_class::get_default_stdout() const{
      return this->default_labels.stdout_name.c_str();
    };
    //-------------------------------------------------------------------//
    const string model_class::get_default_stdmsg() const{
      return this->default_labels.stdmsg_name.c_str();
    };
    //-------------------------------------------------------------------//
    void model_class::set_default_stdin(const char * reassign){
      this->default_labels.stdin_name = reassign;
    };
    //-------------------------------------------------------------------//
    void model_class::set_default_stdmsg(const char * reassign){
      this->default_labels.stdmsg_name = reassign;
    };
    //-------------------------------------------------------------------//
    void model_class::set_default_stdout(const char * reassign){
      this->default_labels.stdout_name = reassign;
    };
    //-------------------------------------------------------------------//
    void model_class::set_default_namelist(const char * reassign){
      this->default_labels.namelist = reassign;
    };
    //-------------------------------------------------------------------//
    void model_class::setdebugonerror(){
      int index = this->get_global_param_index_from_name(gdb_str_DEBUGONERROR);
      this->GLOBAL_PARAMETERS[index]->setValue(1);
    };
    //-------------------------------------------------------------------//
    bool model_class::debugonerror() const{
      return this->get_global_param_by_name(gdb_str_DEBUGONERROR).getiValue() != 0;
    };
    //-------------------------------------------------------------------//
    bool model_class::iscompiled() const{
      return this->model_state.compiled;
    }
    //-------------------------------------------------------------------//
    void model_class::clearspreprocdata(){
      this->spre_data.clear();
    }
    //-------------------------------------------------------------------//
    const spreprocdata & model_class::getspreprocData() const{
      return this->spre_data;
    }
    //-------------------------------------------------------------------//
    model_class::~model_class()
    {
      for(unsigned int i = 0; i < this->GLOBAL_PARAMETERS.size(); i++){
        delete GLOBAL_PARAMETERS[i];
      }
      for(unsigned int i = 0; i < this->GLOBAL_PROFILES.size(); i++){
        delete GLOBAL_PROFILES[i];
      }
      this->MSGOUTs[0].flush();
      this->STDOUTs[0].flush();
      for(size_t iter = 0; iter < this->mdiags.size(); iter++){
        delete this->mdiags[iter];
      }
      mdiags.clear();
    };
    //-----------------------------------------------------
    void model_class::assigndiagfunctions(ParserX & p){  
      assigndiagfuncs(p, true);
      // now go through and assign custom functions
      vector< model_diagnostics_block * > funclist;
      this->getfunctiondiaglist( funclist );
      for(size_t i = 0; i < funclist.size(); i++){
        model_diagnostics_block * dblock = funclist[i];
        p.DefineFun( new user_defined_function( dblock, dblock->GetReturnType(), dblock->GetCallStackTypes()) );
      }
    };
    //-----------------------------------------------------
    void model_class::setbreakpoint(){
      this->mdebugger.setbreakpoint();
    }
    //-----------------------------------------------------
    bool model_class::isbreakpointon(){
      return this->mdebugger.isbreakpointon();
    }
    //-----------------------------------------------------
    void model_class::removebreakpoint(){
      this->mdebugger.removebreakpoint();
    }
    //-----------------------------------------------------
    // Function below sets parameter defaults as well as
    // defines the names and values of parameters needed
    // for the model (also for profiles)
    //-------------------------------------------------------------------//
    void model_class::setAfieldindices(const vector<rank3> & vec){
//#ifdef DEBUG
      if(this->Afieldindices.size() > 0){
        fatal_stop("A field index map has already been assigned to the model.");
      }
//#endif
      this->Afieldindices = vec;
    };
    //-------------------------------------------------------------------//
    const vector<rank3> & model_class::getAfieldindices() const{
      return this->Afieldindices;
    }
    //-------------------------------------------------------------------//
    const ssystem & model_class::getSsystem() const {
      return this->smodule;
    }
 
    //-------------------------------------------------------------------//
    ssystem & model_class::getSsystem(){
      return this->smodule;
    }
    //-------------------------------------------------------------------//
    const psystem & model_class::getPsystem() const{
      return this->pmodule;
    }
    //-------------------------------------------------------------------//
    psystem & model_class::getPsystem(){
      return this->pmodule;
    }
    //-------------------------------------------------------------------//
    qstream & model_class::getSTDMSG() {
#ifdef SERIALQGCM
      return this->MSGOUTs[0];
#else
      if( ! omp_in_parallel() ){ // are we inside a parallel block?
        return this->getmasterSTDMSG();
      }else{ 
        int thread_id = omp_get_thread_num();
        if(thread_id < static_cast<int>(this->out_handle_size)){
          return this->MSGOUTs[thread_id + 1];
        }else{
          fatal_stop("Thread attempted to write msg to handle that was not allocated.\n");
        }
      }
#endif
    };
    //-------------------------------------------------------------------//
    const qstream & model_class::getSTDMSG() const {
#ifdef SERIALQGCM
      return this->MSGOUTs[0];
#else
      if( ! omp_in_parallel() ){ // are we inside a parallel block?
        return this->getmasterSTDMSG();
      }else{ 
        int thread_id = omp_get_thread_num();
        if(thread_id < static_cast<int>(this->out_handle_size)){
          return this->MSGOUTs[thread_id + 1];
        }else{
          fatal_stop("Thread attempted to write msg to handle that was not allocated.\n");
        }
      }
#endif
    };
    //-------------------------------------------------------------------//
    qstream & model_class::getmasterSTDOUT() {
      return this->STDOUTs[0];
    }
    //-------------------------------------------------------------------//
    const qstream & model_class::getmasterSTDOUT() const{
      return this->STDOUTs[0];
    }
    //-------------------------------------------------------------------//
    qstream & model_class::getmasterSTDMSG() {
      return this->MSGOUTs[0];
    }
    //-------------------------------------------------------------------//
    const qstream & model_class::getmasterSTDMSG() const{
      return this->MSGOUTs[0];
    }
    //-------------------------------------------------------------------//
    qstream & model_class::getSTDIN() {
      if(! this->STDIN.reopen() ){ // reopening does not work on internal bufr types
        // so the action will be ignored if this points the stdin
        fatal_stop("Could not reopen stdin for reading.");
      }
      return this->STDIN;
    };
    //-------------------------------------------------------------------//
    bool model_class::setSTDMSG(const char * newfile){
      if(newfile==NULL) return false;
      qstream & stdmsg = this->getSTDMSG();
      stdmsg.close();
      if( ! stdmsg.open(newfile, "w") ){
        return false;
      }
      return true;
    };
    //-------------------------------------------------------------------//
    void model_class::closeDBGIO(){
      this->mdebugger.closeDBGIO();
    };
    //-------------------------------------------------------------------//
    bool model_class::setDBGIN(const char * newfile){
      return this->mdebugger.setDBGIN(newfile);
    };
    //-------------------------------------------------------------------//
    bool model_class::setDBGOUT(const char * newfile){
      return this->mdebugger.setDBGOUT(newfile);
    };
    //-------------------------------------------------------------------//
    bool model_class::setSTDOUT(const char * newfile){
      if(newfile==NULL) return false;
      qstream & stdout = this->getSTDOUT();
      stdout.close();
      if( ! stdout.open( newfile, "w" ) ){
        return false;
      }
      return true;
    };
    //-------------------------------------------------------------------//
    bool model_class::readfifo() {
      const char * fifoname = this->get_global_param_by_name(gdb_str_MODEL_FIFO).getsValue();
      vector<string> input_lines;
      int status = read_fifo(fifoname, input_lines);
      if(status != 0){
        this->printfmsg("Unable to read input from fifo.\n");
        return false;
      };
      if(input_lines.size() == 0) return true;
      // split the input into a vector of lines
      ParserX P(pckNONE);
      assigndiagfunctions(P);
      assignparserdata(P, "v", false);
      for(size_t i = 0; i < input_lines.size(); i++){
        const string & input = input_lines[i];
        if(input.size()>0){
          this->printfmsg("Attempting to evaluate expression '%s'.\n", input.c_str());
          token_list expr;
          extract_words(input, expr);
          P.SetExpr(expr);
          try {
            const string sval = cast_to_string(P.Eval());
            this->printfmsg("Evaluation of fifo expression returned a value of '%f'.\n", sval.c_str());
          }catch(const ParserError & e){
            this->printfmsg("Evaluation of fifo expression returned error message: '%s'.\n", e.what());
            return false;
          }
        }
      }
      return true;
    };
    //-------------------------------------------------------------------//
    const qstream & model_class::getDBGOUT() const{
      return this->mdebugger.getDBGOUT();
      /*if(! this->DBGIO.is_closed()){
        return this->DBGIO;
        }
        return this->getSTDOUT(); */
    };
    //-------------------------------------------------------------------//
    bool model_class::isdebuginputvalid() const {
      return this->mdebugger.isdebugIOvalid();
      /*if( ! this->DBGIO.is_closed() ) return true;
        if( ! this->STDIN.is_closed() ){
        if( STDIN.is_interactive_input() ) return true;
        }
        return false;*/
    };
    //-------------------------------------------------------------------//
    qstream & model_class::getDBGIN() {
      return this->mdebugger.getDBGIN();
      /* if(! this->DBGIO.is_closed()){
         return this->DBGIO;
         }
         return this->getSTDIN(); */
    };
    //-------------------------------------------------------------------//
    const qstream & model_class::getSTDOUT() const{
#ifdef SERIALQGCM
      return this->STDOUTs[0];
#else
      if( ! omp_in_parallel() ){ // are we inside a parallel block?
        return this->getmasterSTDOUT();
      }else{
        int thread_id = omp_get_thread_num();
        if(thread_id < static_cast<int>(this->out_handle_size)){
          return this->STDOUTs[thread_id + 1];
        }else{
          fatal_stop( "Thread attempted to write msg to handle that was not allocated.");
          return stderr;
        }
      }
#endif
    };
    //-------------------------------------------------------------------//
    qstream & model_class::getSTDOUT() {
#ifdef SERIALQGCM
      return this->STDOUTs[0];
#else
      if( ! omp_in_parallel() ){ // are we inside a parallel block?
        return this->getmasterSTDOUT();
      }else{
        int thread_id = omp_get_thread_num();
        if(thread_id < static_cast<int>(this->out_handle_size)){
          return this->STDOUTs[thread_id + 1];
        }else{
          fatal_stop( "Thread attempted to write msg to handle that was not allocated.");
        }
      }
#endif
    };
    //-------------------------------------------------------------------//
    // add macro list function takes a comma separated list of 
    // macro token names that can have definitions and defines them in
    // the model.
    //-------------------------------------------------------------------//
    bool model_class::addmacrolist(const string & defineexpr){
      return this->preprocdata.addmacrolist(defineexpr);
    };
    //-------------------------------------------------------------------//
    const map<string, token_list > &  model_class::getmacros() const 
    {
      return this->preprocdata.getmacrolist();
    }
    //-------------------------------------------------------------------//
    // Macro functions necessary for macro substitution
    //-------------------------------------------------------------------//


    bool model_class::definemacro(const string & name, const token_list & value){
      return this->preprocdata.definemacro(name, value);
    }

    bool model_class::definemacro(const string & name, const string & value){
      return this->preprocdata.definemacro(name, value);
    }

    bool model_class::undefinemacro(const string & name){
      return this->preprocdata.undefinemacro(name);
    };
    //-------------------------------------------------------------------//
    bool model_class::issymboldefined(const char * symbol) const{
      string name(symbol);
      for(unsigned int i = 0; i < pmodule.getmeanfieldsize(); i++){
        if(pmodule.getmeanfieldname(i) == name) return true;
      }

      for(unsigned int i = 0; i < pmodule.getfunctionssize(); i++){
        if(pmodule.getfunctionname(i) == name) return true;
      }

      for(unsigned int i = 0; i < pmodule.getsmeanfieldsize(); i++){
        if(pmodule.getsmeanfieldname(i) == name) return true;
      }

      if( pmodule.doesvariableexist( name ) ) return true;

      if( this->does_prof_exist(symbol) ) return true;

      if( this->does_param_exist(symbol) ) return true;

      if( this->preprocdata.doesmacroexist( symbol ) ) return true;

      return false;
    }
    //-------------------------------------------------------------------//
    void model_class::addparameter(parameter * newparam){
      this->GLOBAL_PARAMETERS.push_back(newparam);
    }
    //-------------------------------------------------------------------//
    void model_class::addprofile(parameter * newprof){
      this->GLOBAL_PROFILES.push_back(newprof);
    }
    //-------------------------------------------------------------------//
    bool model_class::does_param_exist( const char * name ) const{
      for(size_t i = 0; i < this->GLOBAL_PARAMETERS.size(); i++){
        if( this->GLOBAL_PARAMETERS[i]->GetName() == string(name) ) return true;
      }
      return false;
    };
    //-------------------------------------------------------------------//
    bool model_class::does_prof_exist(const char * name) const{
      for(unsigned int i = 0; i < this->GLOBAL_PROFILES.size(); i++){
        if( this->GLOBAL_PROFILES[i]->GetName() == string(name) ) return true;
      }
      return false;
    }
    //-------------------------------------------------------------------//
    int model_class::get_global_param_index_from_name(const char * name) const{
      // this function returns -1 if not found!
      for(unsigned int i = 0; i < this->GLOBAL_PARAMETERS.size(); i++){
        if(strcmp(name, this->GLOBAL_PARAMETERS[i]->GetName()) == 0){
          return i;
        }
      }
      return -1;
    }
    //-------------------------------------------------------------------//
    int model_class::flushmessages() const{
      // flush all output (sync if operating system supports it)
      for(size_t i = 0; i < out_handle_size; i++){
        MSGOUTs[i].flush();
      }
      return 0;
    };
    //-------------------------------------------------------------------//
    int model_class::flushstdout() const{
      // flush all output (sync if operating system supports it)
      for(size_t i = 0; i < out_handle_size; i++){
        STDOUTs[i].flush(); 
      }
      return 0;
    };
    //-------------------------------------------------------------------//
    int model_class::printfmsg(const char * format,...) const{
      va_list(arg);
      va_start(arg, format);
      int rvalue = this->getSTDMSG().vprintf(format, arg);
      return rvalue;
    };
    //-------------------------------------------------------------------//
    int model_class::printf(const char * format, ...) const{
      va_list(arg);
      va_start(arg, format);
      int num = this->getSTDOUT().vprintf(format, arg);
      return num;
    };
    //-------------------------------------------------------------------//
    const parameter& model_class::get_global_param_by_name(const char * name) const{
      for(unsigned int i = 0; i < this->GLOBAL_PARAMETERS.size(); i++){
        if(strcmp(name, (this->GLOBAL_PARAMETERS[i])->GetName()) == 0){
          return *(this->GLOBAL_PARAMETERS[i]);
        }
      }
      string msg;
      msg = "model_class::get_global_param_by_name Error: Could not find parameter: '" + string(name) + "'";
      cgcmexception
        error(msg.c_str(), -4);
      throw error;
      return *(GLOBAL_PARAMETERS[0]);
    }
    //-------------------------------------------------------------------//
    parameter & model_class::get_global_param_by_index(const size_t & i){
      if( i < this->GLOBAL_PARAMETERS.size()){
        return *(this->GLOBAL_PARAMETERS[i]);
      }
      cgcmexception
        error("model_class::get_global_param_by_index Error: Index out of bounds.", -4);
      throw error;
      return (*GLOBAL_PARAMETERS[0]);
    }  
    //-------------------------------------------------------------------//
    const parameter & model_class::get_global_param_by_index(const size_t & i) const{
      if( i < this->GLOBAL_PARAMETERS.size()){
        return *(this->GLOBAL_PARAMETERS[i]);
      }
      cgcmexception
        error("model_class::get_global_param_by_index Error: Index out of bounds.", -4);
      throw error;
      return (*GLOBAL_PARAMETERS[0]);
    }  
    //-------------------------------------------------------------------//
    const parameter & model_class::get_global_prof_by_index(const size_t & i) const{
      if(i < this->GLOBAL_PROFILES.size()){
        return *(this->GLOBAL_PROFILES[i]);
      }
      cgcmexception
        error("model_class::get_global_prof_by_index Error: Index out of bounds.", -4);
      throw error;
      return (*GLOBAL_PROFILES[0]);
    }
    //-------------------------------------------------------------------//
    parameter & model_class::get_global_prof_by_index(const size_t & i){
      if(i < this->GLOBAL_PROFILES.size()){
        return *(this->GLOBAL_PROFILES[i]);
      }
      cgcmexception
        error("model_class::get_global_prof_by_index Error: Index out of bounds.", -4);
      throw error;
      return (*GLOBAL_PROFILES[0]);
    }
    //-------------------------------------------------------------------//
    int model_class::get_global_prof_index_from_name(const char * name) const{
      // this function returns -1 if not found!
      for(unsigned int i = 0; i < this->GLOBAL_PROFILES.size(); i++){
        if(strcmp(name, this->GLOBAL_PROFILES[i]->GetName()) == 0){
          return i;
        }
      }
      return -1;
    }
    //-------------------------------------------------------------------//
    const parameter & model_class::get_global_prof_by_name( const char * name) const{
      for(unsigned int i = 0; i < this->GLOBAL_PROFILES.size(); i++){
        if(strcmp(name, this->GLOBAL_PROFILES[i]->GetName()) == 0){
          return *(this->GLOBAL_PROFILES[i]);
        }
      }
      // include profile name
      string sname(name);
      string speriod(".");
      string errostr = "model_class::get_global_prof_by_name Error: Could not find  requested profile: " + sname + speriod;
      cgcmexception
        error(errostr.c_str(), -4);
      throw error;
      return (*GLOBAL_PROFILES[0]);
    }
    //-------------------------------------------------------------------//
    const string & model_class::getnamelistname(const size_t & index) const{
      const vector<namelistdata*> & namelists = this->preprocdata.namelists;
      if(index >=namelists.size()){
        fatal_stop("Error: Attempted to access element of this->namelists out of bounds.");
      }
      return namelists[index]->filename;
    };
    //-------------------------------------------------------------------//
    INDEX model_class::get_atmosphere_layer_1_field_number() const{
      return this->atmosphere_layer_1;
    };
    //-------------------------------------------------------------------//
    INDEX model_class::get_atmosphere_layer_2_field_number() const{
      return this->atmosphere_layer_2;
    };
    //-------------------------------------------------------------------//
    INDEX model_class::get_ocean_layer_1_field_number() const{
      return this->ocean_layer_1;
    };
    //-------------------------------------------------------------------//
    INDEX model_class::get_ocean_layer_2_field_number() const{
      return this->ocean_layer_2;
    };
    //-------------------------------------------------------------------//
    INDEX model_class::get_mixed_layer_field_number() const{
      return this->mixed_layer;
    };
    //-------------------------------------------------------------------//
    INDEX model_class::get_topo_layer_field_number() const{
      return this->topo_layer;
    };
    //-------------------------------------------------------------------//
    void model_class::removeallbreakpoints(){
      this->mdebugger.removeallbreakpoints();
    };
    //-------------------------------------------------------------------//
    void model_class::getfunctiondiaglist(vector< model_diagnostics_block * > & funclist ){
      for(size_t iter = 0; iter < this->mdiags.size(); iter++){
        if(mdiags[iter]->isblockfunction()){
          funclist.push_back(  mdiags[iter]  );
        }
      }
    };
    //-------------------------------------------------------------------//
    void model_class::getdiaglist(vector< const model_diagnostics_block * > & funclist ) const{
      funclist.resize( mdiags.size(), NULL );
      for(size_t iter = 0; iter < this->mdiags.size(); iter++){
        funclist[iter] =   mdiags[iter];
      }
    }
    //-------------------------------------------------------------------//
    void model_class::getfunctiondiaglist(vector< const model_diagnostics_block * > & funclist ) const{
      for(size_t iter = 0; iter < this->mdiags.size(); iter++){
        if(mdiags[iter]->isblockfunction()){
          funclist.push_back(  mdiags[iter]  );
        }
      }
    };
    //-------------------------------------------------------------------//
    const namelistdata & model_class::getnamelistbyid( size_t i ) const {
	rassert( preprocdata.namelists.size() > i );
    	const namelistdata & ndata = *(this->preprocdata.namelists[i]);
	return ndata; 
    }
    //-------------------------------------------------------------------//
    void model_class::readstring(const char * prompt, string & input) {
      if(this->model_state.io_initialized){
        terminal_readline(input, prompt, this->getSTDIN(), this->getSTDOUT() );
      }else{
        terminal_readline(input, prompt, stdin, stdout);
      };
    };
    //-------------------------------------------------------------------//
    bool model_class::runstatusblocks(bool debug){
      for(size_t i = 0; i < mdiags.size(); i++){
        if(! mdiags[i]->runstatusblock(debug) ) return false;
      }
      return true;
    };
    //-------------------------------------------------------------------//
    bool model_class::rundiagnosticblocks(bool start, bool end){
#ifdef DEBUG
      assert(!(start == true && end == true));
#endif
      if(start){
        for(size_t i = 0; i < mdiags.size(); i++){ if(! mdiags[i]->runstartblock()) return false;}
      }else if(end){ 
        for(size_t i = 0; i < mdiags.size(); i++){ if(! mdiags[i]->runendblock()) return false;}
      }else{
        for(size_t i = 0; i < mdiags.size(); i++){ if(! mdiags[i]->runrepeatblock()) return false;}
      }
      return true;
    }
    //--------------------------------------------------------------------------------------//
    // function below automatically reads a macro substitution and returns
    // a processed string. These macros are either terminated by an end of line
    // or by a semicolon

    //-------------------------------------------------------------------//
    // Function beneath is a parser auxilary function that loads the constant
    // value of all previously defined data into th eparser as well as special
    // tokens such as PI, etc ...
    //-------------------------------------------------------------------//
    void model_class::assignparserparams(ParserX &P) const
    {
      // First assign all auxilary functions used for both namelist processing
      // and the evaluator
      const model_class & mo = *(this);
      // This function stores all parameters as CONSTANTS in the parser
      const size_t psize = mo.get_num_global_params();
      for(size_t i = 0; i < psize; i++){
        const parameter & param = mo.get_global_param_by_index(i);
        if(param.getType() == parameter::Real){
          Value val( static_cast<float_type>( param.getfValue() ) );
          P.DefineConst( param.GetName(), val );
        }else if(param.getType() == parameter::Integer){
          Value val(static_cast<int_type>( param.getiValue() ));
          P.DefineConst( param.GetName(), val );
        }else{ // type is a string
          Value val( param.getsValue() ); 
          P.DefineConst( param.GetName(), val );
        }
      }
    }
    // ---------------------------------------------------------------- //
    void model_class::assignparserprofs(ParserX &P) const
    {
      const model_class & mo = *this;
      for(size_t i = 0; i < mo.get_num_global_profs(); i++){
        const parameter & prof = mo.get_global_prof_by_index(i);
        if( prof.getType() == parameter::Real ){
          Value arg( prof.GetSize(), static_cast<float_type>(0.0) );
          for(size_t j = 0; j < prof.GetSize(); j++) arg[j] = prof.getData()[j];
          P.DefineConst( prof.GetName(), arg );
        }else if( prof.getType() == parameter::Integer ){
          Value arg( prof.GetSize(), static_cast<int_type>(0) );
          for(size_t j = 0; j < prof.GetSize(); j++) arg[j] = prof.getiData()[j];
          P.DefineConst( prof.GetName(), arg );
        }else{ // it's a string
          vector<Value> args(prof.GetSize(), Value(""));
          Value arg(args);
          for(size_t j = 0; j < prof.GetSize(); j++) arg[j] = prof.getsData()[j];
          P.DefineConst( prof.GetName(), arg );
        }
      }
    }
    // ---------------------------------------------------------------- //
    void model_class::assignparserdata(ParserX &P, const char * types, bool optimized, bool diag, int pindex)
    {
      parser_define_basics(P, this, types, optimized, diag, pindex);
      this->assignparserparams(P);
      this->assignparserprofs(P);
      pmodule.assignparserdata(P); // assign P system stuff
    }
    //-------------------------------------------------------------------//
    void model_class::compile(){ // const char * namelist){
      this->preprocess();
      if( ! this->model_state.compiled ){
        compiler::compile( this );
        local::checkprofiles( *this );
        local::checkparams( *this );
        this->model_state.compiled = true;
      }
    };
    //-------------------------------------------------------------------//
    int model_class::assignparam(int pindex, const token_list & valExpr)
    {  // pindex is the index of the parameter that is to be assigned
      parameter & param = this->get_global_param_by_index(pindex);
      if(param.IsConstant()){
        fprintf(stderr, "Parameter %s is a constant and cannot be assigned.\n", param.GetName());
        return -2;
      }
      if(valExpr.empty()){
        fprintf(stderr, "Statement expected after assignment operator.\n");
        return -3;
      }
      parameter::type ptype = param.getType();
      // we expect mistakes, so strict type checking!!
      ParserX P(pckNONE);
      this->assignparserdata( P, "ifbsa", false, false );
      P.SetExpr(valExpr);
      try{
        const IValue & val = P.Eval();
        if(ptype == parameter::Integer){ 
          if(!val.IsInteger() && ( val.GetType() != 'b') ){ 
            fprintf(stderr, "Parameter %s must be assigned an Integer.\n", param.GetName());
            return -2;
          }
          //param.setValue(static_cast<int>(val.GetInteger())); return 0;
          param.setValue( cast_to_int( val ) ); return 0;
        }else if(ptype == parameter::Real){
          if( ( ! val.IsScalar() || val.IsString() || val.IsComplex()) ){
            fprintf(stderr, "Parameter %s must be assigned a Real.\n", param.GetName());
            return -2;	
          }
          param.setValue( static_cast<REAL>( cast_to_float( val ) ) ); return 0;
          // param.setValue( static_cast<REAL>(val.GetFloat()) ); return 0;
        }else if(ptype == parameter::String){
          if( ! val.IsString() ){
            fprintf(stderr, "Parameter %s must be assigned a String.\n", param.GetName()); return -2;
          }
          param.setValue( val.GetString().c_str() );
          return 0;
        }
      }catch(const ParserError & e){
        fprintf(stderr, "Error in expression: %s\n", e.GetContext().Expr.c_str());
        fprintf(stderr, "Error: %s\n", e.what());
        return -2;
      }
      return 0;
    }
    //-------------------------------------------------------------------//
    int assignprofaux(model_class & mo, parameter & prof, const vector<token_list> & vals, const parameter::type ptype){
      const int new_size = vals.size();
      prof.resize( new_size );
      for(size_t i = 0; i < vals.size(); i++){
        ParserX P(pckNONE);
        mo.assignparserdata(P, "ifsa", false, false, i);
        P.DefineConst("this", static_cast<int_type>( i+1 ) );
        P.SetExpr(vals[i]);
        try{
          const IValue * base_val = &(P.Eval());
          if( base_val->GetType() == 'a' ){ // pull out the ith corresponding element 
            const array_type & array = base_val->GetArray();
            if( array.size() > i ){
              base_val = const_cast<Value*>(&(array[i]))->AsIValue();
            } 
          }
          const IValue & val = *base_val;
          //delete [] var_storage; var_storage = NULL;
          // type enforce integer!!
          switch( ptype ){
            case parameter::Integer : 
              if( ! does_value_conform_to_type( val, 'i' ) ) mpError("Cannot convert '%s' to type  '%s'.", convert_expr_to_string(vals[i]).c_str(), get_type_id('i'));
              prof.setValue( static_cast<int>( val.GetInteger() ), i ); break;
            case parameter::Real:
              if( ! does_value_conform_to_type( val, 'f' ) ) mpError("Cannot convert '%s' to type  '%s'.", convert_expr_to_string(vals[i]).c_str(), get_type_id('f'));
              prof.setValue( static_cast<REAL>( val.GetFloat() ), i ); break;
            case parameter::String:
              if( ! does_value_conform_to_type( val, 's' ) ) mpError("Cannot convert '%s' to type  '%s'.", convert_expr_to_string(vals[i]).c_str(), get_type_id('s'));
              prof.setValue( val.GetString().c_str(), i ); break;
            default: fatal_stop("Unrecognized type.");
          };
        }catch(const ParserError & e){
          fprintf(stderr, "Error: %s\n", e.what());
          fprintf(stderr, "Error in expression: %s\n", e.GetContext().Expr.c_str());
          fprintf(stderr, "Error in assignment of element #%i of profile '%s'.\n", static_cast<int>(i)+1, prof.GetName());
          return -2;
        }
      }
      return 0;
    };
    //-------------------------------------------------------------------//
    int assignprofaux_array(model_class & mo, parameter & prof, const vector<token_list> & vals, const parameter::type ptype){
      const int new_size = vals.size();
      try{
        if( new_size != 1) mpError("Right hand side of expression must be an array or list.");
        ParserX P(pckNONE);
        mo.assignparserdata(P, "ifsa", false, false);
        P.SetExpr(vals[0]);
        const IValue & val = P.Eval();
        if( !does_value_conform_to_type( val, 'a' ) ) mpError("Right hand side of expression must be an array or list.");
        const array_type & array = val.GetArray();
        if( array.size() <= 0 ) mpError("Right hand side of expression must be an array or list.");
        prof.resize( array.size() );
        for(size_t i = 0; i < array.size(); i++){ // check types and assign array elements
          const IValue & val = array[i]; // shadow the old val, but no longer needed (lazy copied code here);
          switch( ptype ){
            case parameter::Integer : 
              if( ! does_value_conform_to_type( val, 'i' ) ) mpError("Cannot convert '%s' to type  '%s'.", convert_expr_to_string(vals[i]).c_str(), get_type_id('i'));
              prof.setValue( static_cast<int>( val.GetInteger() ), i ); break;
            case parameter::Real:
              if( ! does_value_conform_to_type( val, 'f' ) ) mpError("Cannot convert '%s' to type  '%s'.", convert_expr_to_string(vals[i]).c_str(), get_type_id('f'));
              prof.setValue( static_cast<REAL>( val.GetFloat() ), i ); break;
            case parameter::String:
              if( ! does_value_conform_to_type( val, 's' ) ) mpError("Cannot convert '%s' to type  '%s'.", convert_expr_to_string(vals[i]).c_str(), get_type_id('s'));
              prof.setValue( val.GetString().c_str(), i ); break;
            default: fatal_stop("Unrecognized type.");
          };
        }
      }catch(const ParserError & e){
        fprintf(stderr, "%s\n", e.what());
        fprintf(stderr, "Error in expression: %s\n", e.GetContext().Expr.c_str());
        fprintf(stderr, "Error in assignment of of profile '%s'.\n", prof.GetName());
        return -1;
      }
      return 0;
    };
    //-------------------------------------------------------------------//
    REAL model_class::computeA(const rank6 & pos) const{
      return Acoeffs.getvalue(pos);
    }

    REAL model_class::computeQ(const rank3 & pos) const{
      return Qcoeffs.getvalue(pos);
    }

    REAL model_class::computeQI(const rank3 & pos) const{
      return QIcoeffs.getvalue(pos);
    }
    //-------------------------------------------------------------------//
    int model_class::assignprof(int pindex, const vector<token_list> & vals, bool incremental, bool array_assignment)
    {
      parameter & prof = this -> get_global_prof_by_index(pindex);
      // Profiles can be of any type now, type checking must be strict!
      const parameter::type ptype = prof.getType();
      int status = 0;
      const vector<token_list> * val_list;
      vector<token_list> mod_val_list;
      if(incremental){
        // need to generate a mod_val_list;   
        val_list = &(mod_val_list);
        // incremental assigments have the form [ init assignment, conditional, step, function iterator];
        // example is [ iter=1, iter < 10, iter=iter+10, sin(iter) ];
        if(vals.size() != 4){
          fprintf(stderr, "Domain assignments require four elements [ incrementing start value, conditional, incrementer, function value ].\n");
          return -4;
        }
        ParserX Pinit(pckNONE ); 
        this->assignparserdata(Pinit, "ibfa", false);
        Value iterator(static_cast<float_type>(0.0)); 
        // get the expression to the left of the equal sign
        const token_list & vartokens = vals[0];
        //preprocessor::tokenizestring(vals[0], vartokens);
        if( vartokens.size() <= 0 ){
          fprintf(stderr, "Domain assignments require a nonempty expression as the first argument to initialize the expression.\n");   return -4;
        }
        const string & itername = vartokens[0]; //.word;
        Pinit.DefineVar( itername, Variable(&iterator) );
        Pinit.SetExpr( vals[0] ); 
        ParserX Pinc = Pinit; Pinc.SetExpr(vals[2]);
        ParserX Pcond = Pinit; Pcond.SetExpr(vals[1]);
        ParserX Pfunc = Pinit; Pfunc.SetExpr(vals[3]);
        try{ 
          for( Pinit.Eval(); cast_to_bool(Pcond.Eval()); Pinc.Eval() ){
            mod_val_list.push_back( token_list(1, cast_to_string( Pfunc.Eval() ) ) ); 
          }
        }catch(const ParserError & e){
          fprintf(stderr, "Error: %s\n", e.what() );
          fprintf(stderr, "Error in expression: %s\n", e.GetContext().Expr.c_str());
          fprintf(stderr, "Error in domain assignment.\n"); return -4;
        }
      }else{
        val_list = &(vals);
      }
      if( array_assignment ){
        status = assignprofaux_array(*this, prof, *val_list, ptype );
      }else{
        status = assignprofaux(*this, prof, *val_list, ptype);
      }
      if( status < 0 ) return status;
      if(strcmp(prof.GetName(), gdb_prof_str_K_DOMAIN) == 0){
        local::mapkcenteredprofile(*this);
        local::mapdkprofile(*this);
      }; 
      return 0;
    }
    //-------------------------------------------------------------------//
    size_t model_class::get_num_global_profs() const {
      return this->GLOBAL_PROFILES.size();
    }
    //-------------------------------------------------------------------//
    size_t model_class::get_num_global_params() const{
      return this->GLOBAL_PARAMETERS.size();
    }
    //-------------------------------------------------------------------//
    bool model_class::nc_save_model_state(const char * templ) const{
      string fnamestr(templ);
      //fnamestr += ".nc"; // add netcdf extension
      const char * fname = fnamestr.c_str();

      nchandler nc_container;
      NcFile::FileMode fmode = NcFile::New;
      double nchandle = 0;

      // the file must not exist, check to make sure it doesn't !!
      nchandle = nc_container.createncfile(fname, NcFile::ReadOnly);
      if(nchandle != 0){
        this->printfmsg("File '%s' already exists. Cannot save model state to that file.\n", fname);
        return false;
      }

      // file does not exist
      nc_container.deletencfile(static_cast<int>(nchandle));
      nchandle = nc_container.createncfile(fname, fmode);
      if(nchandle == 0){ // error occured
        this->printfmsg("File '%s' could not be created.\n", fname);
        return false;
      }
      NcFile * ncfile = nc_container.getncfile(static_cast<int>(nchandle));
      // now we need to go through and save all parameters, profiles
      if(!local::save_model_parameters(*this, ncfile)) return false;
      if(!local::save_model_profiles(*this, ncfile)) return false;
      if(!local::save_model_script_info(*this, ncfile, templ)) return false;

      // add the relevant dimensions
      int tindex = this->pmodule.get_tindex();
      ncfile->add_dim("time", tindex+1);

      // number of fields
      int numfields = this->get_global_param_by_name(gdb_str_NUMFIELDS).getiValue(); 
      ncfile->add_dim("field", numfields);

      // number of fourier angular modes
      const vector<int> & nkdomain = ( this->get_global_prof_by_name(gdb_prof_str_NK_DOMAIN).getiData() );
      ncfile->add_dim("nkindex", static_cast<int>( nkdomain.size() ) );

      // number of kcentered values
      const vector<REAL> & kdomain_centered = (this->get_global_prof_by_name(gdb_prof_str_K_DOMAIN_CENTERED).getData());
      ncfile->add_dim("kindex", static_cast<int>(kdomain_centered.size()));

      if(!pmodule.save_model_state(ncfile)){
        this->printfmsg("Error saving the state of the P-system.\n");
        return false;
      }

      if(!smodule.save_model_state(ncfile)){
        this->printfmsg("Error saving the state of the S-system.\n");
        return false;
      }

      return true;
    };
    //-------------------------------------------------------------------//
    void model_class::nc_restore_model_state_from_file(const char * fname){

      this->init_state(); // initialize to highest level
      nchandler nc_container;
      fprintf(stderr, "========================================================\n");
      fprintf(stderr, "Attempting to restore model state from file '%s'.\n", fname);
      fprintf(stderr, "--------------------------------------------------------\n");
      // the model needs to be at the proper state (compiled, and initialized, all subsystems ! );

      double nchandle = 0;
      // the file must exist, check to make sure it does 
      nchandle = nc_container.createncfile(fname, NcFile::ReadOnly);
      if(nchandle == 0){
        fprintf(stderr, "File '%s' does not exist. Cannot load model state from file.\n", fname);
        throw cgcmexception("Failed to load model state from file.", -1);
      }
      NcFile * ncfile = nc_container.getncfile(static_cast<int>(nchandle));

      if( ! local::restore_check_param_profiles(*this, ncfile) ){
        throw cgcmexception("Failed to load model state from file.", -1);
      }

      if( ! local::restore_check_script_functions(*this, ncfile) ){
        throw cgcmexception("Failed to load model state from file.", -1);
      }

      if( ! this->smodule.restore_model_state(ncfile) ){
        throw cgcmexception("Failed to load model state from file.", -1);
      }

      if( ! this->pmodule.restore_model_state(ncfile) ){
        throw cgcmexception("Failed to load model state from file.", -1);
      }

    };
    //-------------------------------------------------------------------//
#endif
